Flow in a Circular Duct

Table of Contents

1. Introduction
2. Defining the Velocity and Pressure
3. Visualizing the Flow

Introduction

In this notebook we will study how a fluid acts when it flows around a bend opposed to moving through a straight section \ Because the fluid is moving around a bend it will be subject to centripetal force, a ficticious force that acts within a rotating reference frame and keeps the particles moving in circular motion \ Therefore, we will need to use Newton's second law of motion $F={\color[rgb]{0.501961,0.250953,0.010028}m}a$ in order to see how the pressure gradient affects the velocity throughout the bend. Newton published this law along with his others in his 1687 book the Principia and is pictured here:

In order to understand the relationship between pressure and velocity we need Bernoulli's equation which was qualitatively postulated by Daniel Bernoulli in 1738 and formally written by Leonhard Euler in 1752 as $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}$ where ${\color[rgb]{0.315209,0.728565,0.037706}p}_0$ is the total or stagnation pressure, ${\color[rgb]{0.315209,0.728565,0.037706}p}$ is the static pressure, and $\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}$ is the dynamic pressure. Bernoulli (left) and Euler (right) are pictured here:

Sources: Isaac Newton, Daniel Bernoulli, Leonhard Euler

The velocity of a flow moving through a bend can be represented in polar coordinates: $$ \large \require{color}\vec{{\color[rgb]{0.059472,0.501943,0.998465}v}} = {\color[rgb]{0.059472,0.501943,0.998465}v}_{r}\hat{e}_{r} + {\color[rgb]{0.059472,0.501943,0.998465}v}_{\theta}\hat{e}_{\theta} $$ Because the flow is steady the streamlines will follow a circular path, therefore the radial component of velocity is zero, and the velocity will point in the angular direction

# import necessary libraries
import plotly.graph_objects as go
import numpy as np

Defining the Velocity and Pressure

# initialize theta values for graph 
t1=np.linspace(-np.pi/8,np.pi/8,50)
t2=np.linspace(-np.pi/2, np.pi/2,50)

# plot duct and force diagram for small area of flow 
fig = go.Figure()
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')

fig.add_scatter(x=1.25*np.cos(t1), y=1.25*np.sin(t1), line_color='black')
fig.add_scatter(x=1.75*np.cos(t1), y=1.75*np.sin(t1), fill='tonexty', line_color='black', fillcolor = 'steelblue')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(np.pi/8), y=np.array([1.25,1.75])*np.sin(np.pi/8), line_color='black', 
                mode='lines')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(-np.pi/8), y=np.array([1.25,1.75])*np.sin(-np.pi/8), line_color='black', 
                mode='lines')

fig.add_annotation(x=1.2, y=0, text=r'$\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}$', ax=0.5, ay=0, xref='x', 
                   yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=1.8, y=0, text=r"$\require{color} {\color[rgb]{0.315209,0.728565,0.037706}p} + \
                   \frac{d {\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r$", ax=3, ay=0, xref='x', yref='y', 
                   axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=0.6, y=0.75, ax=-0.05, ay=-0.09, xref='x', yref='y', axref='x', ayref='y',arrowhead=5)
fig.add_scatter(x=[0.25], y=[0.5], mode='text', text=r'$r_1$')
fig.add_annotation(x=0.6, y=-1.9, ax=-0.05, ay=0.09, xref='x', yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_scatter(x=[0.25], y=[-0.5], mode='text', text=r'$r_2$')

fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

To find how the velocity changes with radius we need to use Newton's second law of motion: $Force = \require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}a$ \ The only force acting on the flow is pressure since the flow is inviscid and therefore does not experience shear force \ The flow changes direction as it moves around the bend meaning there must be a radial force causing the flow to move in a circular path around the origin \ This force is centripetal and is given by the acceleration $a_{c} = -\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} \hat{r}$ \ To find the mass we can use the fact that mass is equivalent to density times volume, when done per unit width this is $\require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}={\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r)$ \ Then we can use the force diagram above and write out Newton's second Law: $$ \large F=\require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}a=-{\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r) \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} = {\color[rgb]{0.315209,0.728565,0.037706}p} \delta l -\left({\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r\right)\delta l $$ This can then be solved for $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}$: $$ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} = {\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} $$ From this we can see that the normal pressure gradient is what causes the flow to move in a circular path around the bend \ To solve for velocity we will need Bernoulli's equation: $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}$ where $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}$ is a constant\ In order to relate this equation to the pressure gradient equation we just solved for we must take the derivative of Bernoulli's equation so that it involves the pressure gradient \ The derivative should be taken with respect to the radius as that is the only variable that changes velocity and pressure, remembering that this is incompressible flow we get: $\require{color}0=\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}+ {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dr}$, plugging in our previous result for $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}$ we have: $\require{color}0={\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}(r) \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{dr}$ \ Next, we can separate the varaibles and integrate both sides: $\require{color}-\frac{dr}{r} = \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{{\color[rgb]{0.059472,0.501943,0.998465}v}(r)} \longrightarrow {\color[rgb]{0.059472,0.501943,0.998465}v}(r)=\frac{C}{r}$ where C is a constant to be determined by initial conditions \ From this we can see that the velocity is only dependent on the radius so velocity will be constant along at each radius, as we expected \ Using the boundary conditions that $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r_1)={\color[rgb]{0.059472,0.501943,0.998465}v}_0$ we can get that $\require{color}C={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}r_{1}$ therefore $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r)={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}\frac{r_{1}}{r}$ where $r_{1}$ is the inner radius of the bend

# define necessary variables to calculate the speed at each radii
v_0 = 100 # initial speed at r_1
r_1 = 1 # inner radius
r_2 = 2 # outer radius

# calculate the speed between r_1 and r_2
r = np.linspace(r_1, r_2, 9)
v = v_0*r_1*(1/r)

To solve for pressure we must use Bernoulli's equation again\ If we assume that at $r_1$ the static pressure is atmospheric $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_a = 101.325 \ kPa$ then we can find the stagnation pressure, this stagnation pressure remains constant throughout the bend \ To do so we use Bernoulli's equation at $r_1$ as we know both the pressure and velocity at that point, using those values we get: $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_0 = {\color[rgb]{0.315209,0.728565,0.037706}p}_{\color[rgb]{0.989013,0.435749,0.811750}a} + \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2$ Now we can solve for the pressure at an radius by using the stagnation pressure we just solved for and the equation for velocity that we previously solved for: $$ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}(r) = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2(r) = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_0^2r_1^2}{r^2} \longrightarrow {\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_a + \frac{1}{2} {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2\left(1-\frac{r_1^2}{r^2}\right) $$

# define variables needed to calculate pressure
rho = 1.293 # denisty of air in kg/m^3
p_a = 101325 # pressure of air (atmospheric pressure) in Pa
p_0 = p_a + (1/2)*rho*v_0**2

# calculate pressure between r_1 and r_2
p = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(r**2)

# define variables needed to make a contour plot
theta = np.linspace(0, np.pi/2, 50)
x=r_2*np.cos(theta)
y=np.linspace(-2,2,50)

Visualizing the Flow

# define pressure and velocity for contour plot
v_z = []
p_z = []
for y_val in y:
    v_y_list = []
    p_y_list = []
    for x_val in x:
        if (x_val**2 + y_val**2) < r_1:
            p_xy = p_a
            v_xy = v_0
        else:
            v_xy = v_0/np.sqrt(x_val**2 + y_val**2)
            p_xy = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(x_val**2 + y_val**2)
        v_y_list.append(v_xy)
        p_y_list.append(p_xy)
    v_z.append(v_y_list)
    p_z.append(p_y_list)

# make contour plots for velocity and pressure  
fig = go.Figure()
fig.add_contour(z=v_z, x=x, y=y, contours_coloring='heatmap', line_width=0, colorbar=dict(title='Velocity (m/s)', 
                titleside='right'), colorscale='PuBuGn', reversescale=True)
fig.add_contour(z=p_z, x=x, y=y, contours_coloring='heatmap', line_width=0, visible=False, colorbar=dict(title=
                'Pressure (Pa)', titleside='right'), colorscale='dense', reversescale=True)

fig.add_shape(type='circle', xref='x', yref='y', x0=-1, y0=-1, x1=1, y1=1, fillcolor='white', line_color='white')
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), line_width=0)
fig.add_scatter(x=3*np.cos(t2), y=3*np.sin(t2), line_color='white', fill='tonexty', fillcolor='white')

fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False, range=[-2,2]), 
                  xaxis=dict(showticklabels=False), plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=800, height=800, 
                  showlegend=False, updatemenus=[dict(active=0, buttons=[dict(label='Velocity', method='update', 
                  args=[{'visible':[True, False, True, True, True]}]), dict(label='Pressure', method='update',
                  args=[{'visible':[False, True, True, True, True]}])], y=1.2, x=0.41)])
fig.show()

As you can see in the contour plots above the pressure and velocity are inversely proportional to eachother which is what we would expect based on Bernoulli's equation\ Additionally, the velocity decreases with increasing radius which we would also expect because centripetal acceleration decreases as the radius of curvative increases

# define lines where we want to visualize the flow and time and distance information for those lines
total_distance = np.pi*r
total_time = total_distance / v
dist = total_time[7] * v

r_lines = np.array([r[1], r[4], r[7]])
r_dist = np.array([dist[1], dist[4], dist[7]])
rad = -np.pi/2 + (r_dist / r_lines)

# create animated plot to show ideallized flow movement 
fig = go.Figure()

frames = [go.Frame(data=[]) for k in range(50)]
for i, r_line in enumerate(r_lines):
    fig.add_scatter(x=r_line*np.cos(t2), y=r_line*np.sin(t2), line_dash='dash', line_color='lightsteelblue')
    x_pos = r_line*np.cos(np.linspace(-np.pi/2, rad[i], 50))
    y_pos = r_line*np.sin(np.linspace(-np.pi/2, rad[i], 50))
    try:
        end = np.where(x_pos < 0)[0][0]
        x_pos[end:len(x_pos)] = 'NaN'
    except IndexError:
        None
    for j, f in enumerate(frames):
        f.data += (go.Scatter(x=[x_pos[j]], y=[y_pos[j]], mode='markers', marker=dict(color='teal', size=10)),)

fig.update(frames=frames)
fig.update_layout(updatemenus = [dict(type='buttons', buttons=[dict(label='play', method='animate', 
                  args=[None, {'frame':{'duration':100}}])])], xaxis=dict(range=[0,2], showticklabels=False), 
                  yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), showlegend=False, 
                  plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=900, height=700)

fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')

fig.show()

In the above animation we can visually see how an idealized fluid would move around a bend, and we can see how the particles at inner radii move faster than those at outer radii